program main
	
	use MI
	
	implicit NONE
	
	integer(4), parameter			:: NINDXESMAX = 2*LMAX+LX*LX
	integer(4), parameter			:: nterm = 50000
	integer(4), parameter			:: nmax = 4 * 131072
	
	integer(4)					:: LC
	
	integer(4)					:: nmc
	real(8)						:: beta, t, s, sq
	integer(4)					:: iunit
	integer(4)					:: H(3, LMAX)
	integer(4)					:: ip
	integer(4)					:: nonzero

	integer(4)					:: ix(NINDXESMAX)
	integer(4)					:: nindxes
	integer(4)					:: inac, kx
	
	real(8)						:: ncdat(nmax)
	
	real(8)						:: sz(LX,LX), szstep(LX,LX)
	integer(4)					:: u1, u2
	
	t = 10.0d0
	
		beta = 1.0d0 / t
		LC = 5
		nindxes = 2*LC+LX*LX
		
		call Init(LC) 
		H(:,1:LC) = 0
		
		call InitAntiFerroSites()
		sites( LX/2 + 1, LX/2 + 1) = sites( LX/2 + 1, LX/2 + 1) * 3
		
		!sz = 0.0d0
		!do inac = 1, nterm
		!	call GlobalStep(LC,H(:,1:LC))
		!	call CalcSz( szstep, LC, H)
		!	sz = sz + szstep
		!end do
		!sz = sz / nterm
		!print *, sz!( LX/2 + 1, LX/2 + 1)
		!stop
		
		do inac = 1, nterm
			call SetRandomSequence(nindxes,ix(1:nindxes))
			do kx = 1, nindxes
				
				if ( ix(kx) <= LX*LX) then
					call GlobalStep(LC,H(:,1:LC))
				else if ( ix(kx) <= LX*LX + LC) then
					call LocalDuostep(beta,LC,H(:,1:LC))
				else
					call LocalUniStep(beta,LC,H(:,1:LC))
				end if
				
			end do
			
			nonzero = count( H(1,1:LC) /= 0)
			if ( nonzero > 0.95d0 * LC) then
				if (LC+5 > LMAX) stop "AZAZA"
				H(:,LC+1:LC+5) = 0
				LC = LC + 5
				call Init(LC)
				nindxes = 2*LC+LX*LX
				print *, inac, LC
			end if
			
		end do
		
		s = 0
		sq = 0

		do nmc = 1, nmax
			
			call SetRandomSequence(NINDXES,ix)
			do kx = 1, NINDXES
				if ( ix(kx) <= LX*LX) then
					call GlobalStep(LC,H(:,1:LC))
				else if ( ix(kx) <= LX*LX + LC) then
					call LocalDuostep(beta,LC,H(:,1:LC))
				else
					call LocalUniStep(beta,LC,H(:,1:LC))
				end if
			end do
			
			nonzero = count( H(1,1:LC) /= 0)
			ip = ip + 1
			s = s + nonzero
			sq = sq + nonzero**2
			
			ncdat(nmc) = s / nmc

		end do
		
	open( newunit = iunit, file = 'nc.bin', access = 'stream', action = 'write', status = 'replace')
	write( iunit) ncdat
	close( iunit)
	stop
end program main